function [G1,CC,impact,fmat,fwt,ywt,gev,eu]=amagensys(g0,g1,cc,psi,pi,uprbnd,condn)
%========================================================================== 
% AMAGENSYS.m 
% This code receives Gensys inputs and computes the solution using the 
% Anderson-Moore algorithm and the Gensys to AMA interface 
%
% function [G1,CC,impact,fmat,fwt,ywt,gev,eu]=amagensys(g0,g1,cc,psi,pi,uprbnd,condn); 

% Inputs: 
% -------
% Write the matrices of the standard Gensys algorithm 
% g0*y(t)=g1*y(t-1)+c+psi*z(t)+pi*eta(t),
% Additional inputs from standard gensys used in SPAmalg.m 
% UPRBND    upper bound allowed for unit roots, similar to DIV in
%           gensys
% CONDN     Zero tolerance used as a condition number test
%           by numeric_shift and reduced_form
% If not defined will be set to their default values 
%
% Outputs
% -------
% Same as Gensys solution 
% y(t)=G1*y(t-1)+C+impact*z(t)+ywt*inv(I-fmat*inv(L))*fwt*z(t+1) 
%
% Note: In gensys GEV=[A B] and compute the roots as B./A 
% AMA returns only a subset of the roots so to make output comparable 
% will report an artificial variable GEV=[Aama Bama] 
% where Bama has in the first K rows the roots returned by AMA 
% and zeros in the next (NY-NK) rows, while A is an NY x 1 column of 1s
%
% Alejandro Justiniano
% Adapted from Gary Anderson's GENSYStoAMA matlab code. 
% Changed sorting of output by AIMCODE values to guarantee output is
% produced
% July 9th 2007 
% =========================================================================
if nargin < 7
    condn=1.e-10;
    if nargin < 6
        uprbnd =1+1.e-6;%allow unit roots
    end
end
ny=size(g0,1); 
% 1. Convert Gensys inputs into AMA
% ----------------------------------
[theHM,theH0,theHP]=convertFromGensysIn(g0,g1,pi);
theH=[theHM,theH0,theHP];
% 2, Call SPAMAlg that implements AMA
% -----------------------------------
[bb,rts,ia,nexact,nnumeric,lgroots,aimcode]=SPAmalg(theH,size(theHM,1),1,1,condn,uprbnd);
if aimcode > 5
    G1=[];CC=[];impact=[];ywt=[];fmat=[];fwt=[];gev=[];
    eu=[0;0]; return;
end
eu=setEu(aimcode);eu=eu(:); 
% 3. If a solution exists, convert output from AMA to Gensys 
% -----------------------------------------------------------
if isequal(eu,[1;1])
    phi=(theH0+theHP*sparse(bb))\eye( size(theH0,1));
    [CC,G1,impact,ywt,fmat,fwt]=convertToGensysOut(bb,phi,(-phi*theHP),cc,g0,g1,psi,size(pi,2));
    gev=[ones(ny,1) [rts;ones(ny-length(rts),1)]]; 
else
    %%disp('No unique stable solution with AMA'); 
    G1=[];CC=[];impact=[];ywt=[];fmat=[];fwt=[];gev=[];   
end
% ======================================
% Solution indicator same as Gensys 
% =======================================
function euVal=setEu(aimCode)
euVal=zeros(1,2);
switch aimCode
    case 1
        euVal(1,1)=1;
        euVal(1,2)=1;
    case 2
        euVal(1,1)=-2;
        euVal(1,2)=-2;
        disp('Coincident Zeros');
    case 3
        euVal(1,1)=0;
        euVal(1,2)=0;
        disp('No solution'); 
    case 4
        euVal(1,1)=1;
        euVal(1,2)=0;
        disp('Indeterminacy'); 
    case 5
        euVal(1,1)=-5;
        euVal(1,2)=-5;
end